##### Climate Policy Attitudes in Consumers #####
### March 2025
### Nathan Mariano 

#### Set-up #####
library(readr)
library(plyr)
library(tidyverse)
library(tidylog)
library(ggplot2)
library(broom) 
library(labelled) 
library(kableExtra)
library(car)
library(stargazer)

setwd("/Users/nmariano/Library/CloudStorage/Dropbox/1 Projects/Consumer Preferences/")

surv <- read_csv('Data/consumers_data.csv')

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  require(dplyr)
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- plyr::ddply(data, groupvars, .drop=.drop,
                       .fun = function(xx, col) {
                         c(N    = length2(xx[[col]], na.rm=na.rm),
                           mean = mean   (xx[[col]], na.rm=na.rm),
                           sd   = sd     (xx[[col]], na.rm=na.rm)
                         )
                       },
                       measurevar
  )
  # Rename the "mean" column    
  datac <- plyr::rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  return(datac)
}

#### Clean and Split Data ####
# Split datasets for RMD and GPV analyses
gpv_data <- 
  surv |>
  filter(gpv_treat == 1 | control_group == 1) |>
  mutate(
    party_id = as.factor(party_id),
    Sex =as.factor(Sex),
    ethnicity = as.factor(ethnicity),
    age_group = as.factor(age_group),
    area_of_res = as.factor(area_of_res),
    
    party_id = relevel(party_id, ref = "Republican"), 
    Sex = relevel(Sex, ref = "Female"), 
    ethnicity = relevel(ethnicity, ref = "White"), 
    age_group = relevel(age_group, ref = "36-50"), 
    area_of_res = relevel(area_of_res, ref = "Suburban (residential area surrounding a large city)")
  )

rm_data <- 
  surv |>
  filter(rm_treat == 1 | control_group == 1) |>
  mutate(
    party_id = as.factor(party_id),
    Sex =as.factor(Sex),
    ethnicity = as.factor(ethnicity),
    age_group = as.factor(age_group),
    area_of_res = as.factor(area_of_res),
    
    party_id = relevel(party_id, ref = "Republican"), 
    Sex = relevel(Sex, ref = "Female"), 
    ethnicity = relevel(ethnicity, ref = "White"), 
    age_group = relevel(age_group, ref = "36-50"), 
    area_of_res = relevel(area_of_res, ref = "Suburban (residential area surrounding a large city)")
  )

#### Correlates of RMD Coefplot ####
demo_rm_strong <- 
  lm(rm_cons_strong ~ 
       party_id +
       Sex + 
       ethnicity + 
       area_of_res + 
       age_group, 
     data = rm_data)

demo_rm_strong <- tidy(demo_rm_strong, conf.int = TRUE)

demo_rm_strong_coef <- 
  demo_rm_strong |>
  filter(term != "(Intercept)") |>
  mutate(term = case_when(
    term == 'party_idDemocrat' ~ 'Democrat', 
    term == 'party_idIndependent' ~ 'Independent',
    term == 'SexMale' ~ 'Male', 
    term == 'ethnicityAsian' ~ 'Asian', 
    term == 'ethnicityBlack' ~ 'Black', 
    term == 'ethnicityOther' ~ 'Other', 
    term == 'area_of_resUrban (apartment building or condo in a city)' ~ 'Urban', 
    term == 'area_of_resSemi-urban (residential part of a large city)' ~ 'Semi-urban', 
    term == 'area_of_resRural' ~ 'Rural', 
    term == 'age_group18-35' ~ '18-35', 
    term == 'age_group51-65' ~ '51-65', 
    term == 'age_group66+' ~ '66+'
  )) |>
  filter(
    term %in% c(
      'Democrat', 'Independent', 'Male', 'Asian', 'Black', 'Other',  
      'Urban', 'Semi-urban', 'Rural', '18-35', '51-65', '66+')
  ) |>
  mutate(term = factor(term, 
                       levels = rev(c(
                         "Democrat", 
                         "Independent", 
                         "Male", 
                         "Black", 
                         "Asian", 
                         "Other",
                         "Urban",
                         "Rural",
                         "Semi-urban",
                         "18-35", 
                         "51-65", 
                         "66+")),
                       ordered = TRUE)) |>
  mutate(
    category = case_when(
      term == 'Democrat' ~ 'Party ID', 
      term == 'Independent' ~ 'Party ID', 
      term == 'Male' ~ 'Gender', 
      term == 'Black' ~ 'Ethnicity', 
      term == 'Asian' ~ 'Ethnicity',
      term == 'Other' ~ 'Ethnicity', 
      term == 'Urban' ~ 'Area of Residence', 
      term == 'Semi-urban' ~ 'Area of Residence', 
      term == 'Rural' ~ 'Area of Residence', 
      term == '18-35' ~ 'Age', 
      term == '51-65' ~ 'Age', 
      term == '66+' ~ 'Age'
    )) |>
  mutate(category = factor(category,
                           levels = rev(c(
                             "Age",
                             "Area of Residence",
                             "Ethnicity",
                             "Gender",
                             "Party ID"
                           )),
                           ordered = TRUE)
  ) 


# Make a plot 
ggplot(demo_rm_strong_coef, aes(x = estimate, y = term)) +
  geom_point(position=position_dodge(.4), size=3) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high),
                width=.2,
                position=position_dodge(.4)) +
  # Add vertical dotted line at x = 0
  geom_vline(aes(xintercept = 0), linetype = "dotted") +
  # Adjust labels and titles
  labs(x = "Correlation to RMD Consumption", 
       y = NULL) +
  theme_bw() +
  facet_grid(category ~ ., scales = "free_y", space = "free") +
  theme(
    strip.text.y = element_text(angle = 0) # Optional: rotate category labels for readability
  )

#### RMD Subgroup Coefficient Plot #### 
support_all <- summarySE(data=rm_data, measurevar=c("policy_ft"), groupvars=c("rm_treat"), na.rm=TRUE)
support_rm <- summarySE(data=rm_data[which(rm_data$rm_cons_strong == 1),], measurevar=c("policy_ft"), groupvars=c("rm_treat"), na.rm=TRUE)
support_nonrm <- summarySE(data=rm_data[which(rm_data$rm_cons_strong == 0),], measurevar=c("policy_ft"), groupvars=c("rm_treat"), na.rm=TRUE)

support_graph <- 
  rbind(support_all, support_rm, support_nonrm) |>
  mutate(
    group = c("All", "All", "RMD Cons.", "RMD Cons.", "Non-RMD Cons.", "Non-RMD Cons."), 
    rm_treat = case_when(
      rm_treat == 0 ~ "Control", 
      rm_treat == 1 ~ "Treatment"
    )
  )
support_graph$group <- c("All", "All", "RMD Cons.", "RMD Cons.", "Non-RMD Cons.", "Non-RMD Cons.")

ggplot(support_graph, aes(x=as.factor(rm_treat), y=policy_ft, color=group)) + 
  geom_point(position=position_dodge(.4), size=3) +
  scale_y_continuous(limits=c(0,100)) +
  scale_color_manual(values=c("#000000", "#B2DF8A", "#FB8072")) +
  ylab("Support for Climate Policy") +
  xlab("Treatment Status") +
  theme_bw() +
  geom_errorbar(aes(ymin=policy_ft-ci, ymax=policy_ft+ci),
                width=.2,
                position=position_dodge(.4)) +
  labs(color="Consumer Subgroup") +  # Update legend label
  theme(
    legend.position = c(0.03, 0.97),  # Top-left position within axis area
    legend.justification = c(0, 1),  # Aligns legend to the top-left
    legend.background = element_blank(), # Removes the solid background
    legend.key = element_blank(),      # Removes the background behind keys
    legend.text = element_text(size=8) # Adjusts legend text size
  )

#ggsave('Results/ERL Version/subgroup_coefplot.pdf', width = 7, height = 6, dpi = 300)

#### RMD Demo Controls Coefficient Plot ####
support_controls <- 
  lm(policy_ft ~ 
       rm_cons_strong*rm_treat +
       party_id +
       Sex + 
       ethnicity + 
       area_of_res + age, 
     data = rm_data)

support_controls <- tidy(support_controls, conf.int = TRUE)

support_controls_coef <- 
  support_controls |>
  filter(term != "(Intercept)") |>
  mutate(effect_type = 'Interaction Effects') |>
  mutate(term = case_when(
    term == 'rm_cons_strong:rm_treat' ~ 'RMD Consumer x Policy',
    term == 'rm_treat' ~ "RMD Policy", 
    term == 'rm_cons_strong' ~ 'RMD Consumer',
    term == 'party_idDemocrat' ~ 'Democrat', 
    term == 'SexMale' ~ 'Male', 
    term == 'ethnicityAsian' ~ 'Asian', 
    term == 'ethnicityBlack' ~ 'Black', 
    term == 'area_of_resRural' ~ 'Rural', 
    term == 'age' ~ 'Age'
  )) |>
  filter(
    term %in% c(
      'RMD Consumer x Policy', 'RMD Policy', 'RMD Consumer', 
      'Democrat', 'Male', 'Asian', 'Black', 'Rural', 'Age')
  ) |>
  mutate(term = factor(term, 
                       levels = rev(c(
                         "RMD Consumer x Policy",
                         "RMD Consumer", 
                         "RMD Policy", 
                         "Democrat", 
                         "Male", 
                         "Black", 
                         "Asian", 
                         "Rural", 
                         "Age")),
                       ordered = TRUE))

# Make a plot 
ggplot(support_controls_coef, aes(x = estimate, y = term)) +
  # Add coefficient points
  geom_point(position=position_dodge(.4), size=3) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high),
                width=.2,
                position=position_dodge(.4)) +
  # Add vertical dotted line at x = 0
  geom_vline(aes(xintercept = 0), linetype = "dotted") +
  # Adjust labels and titles
  labs(x = "Effect on Support for Climate Policy", 
       y = NULL) +
  theme_bw()


#### RMD Democrat Interaction Coefficient Plot ####
dem_inter <- 
  lm(policy_ft ~ 
       rm_cons_strong*rm_treat + rm_treat*dem, 
     data = surv)

dem_inter_coef <- tidy(trip_mod, conf.int = TRUE)

dem_inter_coef <- 
  dem_inter_coef |>
  filter(term != "(Intercept)") |>
  mutate(term = case_when(
    term == 'rm_cons_strong:rm_treat' ~ "RMD Consumer x Policy",
    term == 'rm_treat' ~ "RMD Policy", 
    term == 'rm_cons_strong' ~ 'RMD Consumer',
    term == 'dem' ~ 'Democrat', 
    #term == 'rm_cons_strong:rm_treat:dem' ~ 'RMD Consumer x Policy x Dem', 
    term == 'rm_treat:dem' ~ 'RMD Policy x Dem', 
  )) |>
  mutate(term = factor(term, 
                       levels = rev(c(
                         "RMD Consumer x Policy", 
                         "RMD Policy x Dem",
                         "RMD Policy", 
                         "RMD Consumer",
                         'Democrat')),
                       ordered = TRUE))

# Make a plot 
ggplot(dem_inter_coef, aes(x = estimate, y = term)) +
  # Add coefficient points
  geom_point(position=position_dodge(.4), size=3) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high),
                width=.2,
                position=position_dodge(.4)) +
  # Add vertical dotted line at x = 0
  geom_vline(aes(xintercept = 0), linetype = "dotted") +
  # Adjust labels and titles
  labs(x = "Effect on Support for Climate Policy", 
       y = NULL) +
  theme_bw()

#### GPV Subgroup Coefficient Plot ####

support_all <- summarySE(data=gpv_data, measurevar=c("policy_ft"), groupvars=c("gpv_treat"), na.rm=TRUE)
support_gpv <- summarySE(data=gpv_data[which(gpv_data$gpv_cons_strong == 1),], measurevar=c("policy_ft"), groupvars=c("gpv_treat"), na.rm=TRUE)
support_nongpv <- summarySE(data=gpv_data[which(gpv_data$gpv_cons_strong == 0),], measurevar=c("policy_ft"), groupvars=c("gpv_treat"), na.rm=TRUE)

support_graph <- 
  rbind(support_all, support_gpv, support_nongpv) |>
  mutate(
    group = c("All", "All", "GPV Cons.", "GPV Cons.", "Non-GPV Cons.", "Non-GPV Cons."), 
    gpv_treat = case_when(
      gpv_treat == 0 ~ "Control", 
      gpv_treat == 1 ~ "Treatment"
    )
  )
support_graph$group <- factor(support_graph$group, levels = c("All", "Non-GPV Cons.", "GPV Cons."))


ggplot(support_graph, aes(x=as.factor(gpv_treat), y=policy_ft, color=group)) + 
  geom_point(position=position_dodge(.4), size=3) +
  scale_y_continuous(limits=c(0,100)) +
  scale_color_manual(values=c("#000000", "#B2DF8A", "#FB8072")) +
  ylab("Support for Climate Policy") +
  xlab("Treatment Status") +
  theme_bw() +
  geom_errorbar(aes(ymin=policy_ft-ci, ymax=policy_ft+ci),
                width=.2,
                position=position_dodge(.4)) +
  labs(color="Consumer Subgroup") +  # Update legend label
  theme(
    legend.position = c(0.03, 0.97),  # Top-left position within axis area
    legend.justification = c(0, 1),  # Aligns legend to the top-left
    legend.background = element_blank(), # Removes the solid background
    legend.key = element_blank(),      # Removes the background behind keys
    legend.text = element_text(size=8) # Adjusts legend text size
  )

#ggsave('Results/ERL Version/subgroup_coefplot_gpv.pdf', width = 7, height = 6, dpi = 300)

#### GPV Demo Controls Coefficient Plot ####
support_controls <- 
  lm(policy_ft ~ 
       gpv_cons_strong*gpv_treat +
       party_id +
       Sex + 
       ethnicity + 
       area_of_res + age, 
     data = gpv_data)

support_controls <- tidy(support_controls, conf.int = TRUE)

support_controls_coef <- 
  support_controls |>
  filter(term != "(Intercept)") |>
  mutate(effect_type = 'Interaction Effects') |>
  mutate(term = case_when(
    term == 'gpv_cons_strong:gpv_treat' ~ 'GPV Consumer x Policy',
    term == 'gpv_treat' ~ "GPV Policy", 
    term == 'gpv_cons_strong' ~ 'GPV Consumer',
    term == 'party_idDemocrat' ~ 'Democrat', 
    term == 'SexMale' ~ 'Male', 
    term == 'ethnicityAsian' ~ 'Asian', 
    term == 'ethnicityBlack' ~ 'Black', 
    term == 'area_of_resRural' ~ 'Rural', 
    term == 'age' ~ 'Age'
  )) |>
  filter(
    term %in% c(
      'GPV Consumer x Policy', 'GPV Policy', 'GPV Consumer', 
      'Democrat', 'Male', 'Asian', 'Black', 'Rural', 'Age')
  ) |>
  mutate(term = factor(term, 
                       levels = rev(c(
                         "GPV Consumer x Policy",
                         "GPV Consumer", 
                         "GPV Policy", 
                         "Democrat", 
                         "Male", 
                         "Black", 
                         "Asian", 
                         "Rural", 
                         "Age")),
                       ordered = TRUE))

# Make a plot 
ggplot(support_controls_coef, aes(x = estimate, y = term)) +
  # Add coefficient points
  geom_point(position=position_dodge(.4), size=3) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high),
                width=.2,
                position=position_dodge(.4)) +
  # Add vertical dotted line at x = 0
  geom_vline(aes(xintercept = 0), linetype = "dotted") +
  # Adjust labels and titles
  labs(x = "Effect on Support for Climate Policy", 
       y = NULL) +
  theme_bw()

# Save interaction plot 
#ggsave('Results/ERL Version/support_controls_coefplot_gpv.pdf', width = 8, height = 5, dpi = 300)

#### Weak RMD Subgroup Coefficient Plot #### 
support_all <- summarySE(data=rm_data, measurevar=c("policy_ft"), groupvars=c("rm_treat"), na.rm=TRUE)
support_rm <- summarySE(data=rm_data[which(rm_data$rm_cons_weak == 1),], measurevar=c("policy_ft"), groupvars=c("rm_treat"), na.rm=TRUE)
support_nonrm <- summarySE(data=rm_data[which(rm_data$rm_cons_weak == 0),], measurevar=c("policy_ft"), groupvars=c("rm_treat"), na.rm=TRUE)

support_graph <- 
  rbind(support_all, support_rm, support_nonrm) |>
  mutate(
    group = c("All", "All", "Weak RMD Cons.", "Weak RMD Cons.", "Non-RMD Cons.", "Non-RMD Cons."), 
    rm_treat = case_when(
      rm_treat == 0 ~ "Control", 
      rm_treat == 1 ~ "Treatment"
    )
  )
support_graph$group <- c("All", "All", "Weak RMD Cons.", "Weak RMD Cons.", "Non-RMD Cons.", "Non-RMD Cons.")

ggplot(support_graph, aes(x=as.factor(rm_treat), y=policy_ft, color=group)) + 
  geom_point(position=position_dodge(.4), size=3) +
  scale_y_continuous(limits=c(0,100)) +
  scale_color_manual(values=c("#000000", "#B2DF8A", "#FB8072")) +
  ylab("Support for Climate Policy") +
  xlab("Treatment Status") +
  theme_bw() +
  geom_errorbar(aes(ymin=policy_ft-ci, ymax=policy_ft+ci),
                width=.2,
                position=position_dodge(.4)) +
  labs(color="Consumer Subgroup") +  # Update legend label
  theme(
    legend.position = c(0.03, 0.97),  # Top-left position within axis area
    legend.justification = c(0, 1),  # Aligns legend to the top-left
    legend.background = element_blank(), # Removes the solid background
    legend.key = element_blank(),      # Removes the background behind keys
    legend.text = element_text(size=8) # Adjusts legend text size
  )

ggsave('Results/ERL Version/subgroup_coefplot_weak.pdf', width = 7, height = 6, dpi = 300)